-
Notifications
You must be signed in to change notification settings - Fork 47
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Phonopy: fixing atoms #2073
Phonopy: fixing atoms #2073
Conversation
Can one of the admins verify this patch? |
This is worth thinking about. Let me think out loud here. In Scenario A, the user wants to relax all the atoms and then, afterwards, do a phonon calculation with a subset of atoms fixed. In this scenario, they should just run a from quacc.recipes.emt.core import relax_job
from quacc.recipes.emt.phonons import phonon_flow
output = relax_job(atoms)
new_atoms = output["atoms"]
new_atoms.set_constraint(...)
output2 = phonon_flow(new_atoms, run_relax=False) In Scenario B, the user wants to have a subset of atoms fixed during both the relax and the phonon flow. In this case, the only need to call So, in both scenarios, this is currently possible. While I agree that adding a dedicated keyword argument is more convenient for the user, the counterpoint to this is that it adds more lines to maintain even though it is not strictly necessary. Right now, virtually every recipe in quacc supports fixing atoms in any way you'd like simply by passing in an
Oops, yeah that was a mistake. I also like the addition you have made regarding symmetrization of the force constants. |
To use ASE's constraints the only way that comes to my mind is to do: fixed = []
for constr in atoms.constraints:
if isinstance(constr, FixAtoms):
fixed.extend(constr.get_indices()) Then remove the |
@tomdemeyere: Thanks. My point was more about the need for these changes at all if one can pass in an if len(fixed) > 0:
fixed = get_phonopy_structure(AseAtomsAdaptor.get_structure(fixed))
fixed = phonopy.structure.cells.get_supercell(fixed, supercell_matrix)
fixed = AseAtomsAdaptor.get_atoms(get_pmg_structure(fixed)) I don't yet understand what the issues are (in The other changes, specifically regarding getting rid of the second |
Apologies for the lack of context. What we have to do to effectively fix some atoms when using phonopy is the following:
ASE's constraints unfortunately are not of help here. This whole trick has to be done manually. This whole process drastically reduces the number of displacement required, for an adsorbate on a slab for example you might end up calculating the displacement of the adsorbate only. As you might guess this is a huge approximation because you assume that modes from the fixed system do not interact with the ones of the non-fixed system. On the previous MR I did a test which showed that for water on Pt(111) the trick does not change frequencies much as the slab and waters are loosely coupled. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@tomdemeyere: Thank you very much for your helpful reply. I had (incorrectly) assumed based on your PR description that this was mainly centered around ease-of-use. However, I now understand that it is actually correcting the behavior for fixed atoms and that simply supplying the Atoms
object with FixAtoms
constraints is not sufficient. I am familiar with the typical assumption about decoupling of vibrational modes between surface and adsorbate but have normally used ASE's Vibrations
module and ASE's thermo modules for this instead of phonon calculations via phonopy. If I recall correctly, ASE's Vibrations
module is smart enough not to displace and vibrate atoms that are fixed.
Anyway, I have left several comments below. Apologies in advance about some of them being a bit pedantic, but it will help increase the clarity for future me.
src/quacc/recipes/common/phonons.py
Outdated
fixed_atoms = np.full(len(phonon.supercell), False) | ||
fixed_atoms = np.append(fixed_atoms, [True] * len(atoms_to_add)) | ||
fixed_atoms = fixed_atoms.astype(bool) | ||
|
||
supercells = [ | ||
phonopy_atoms_to_ase_atoms(s) for s in phonon.supercells_with_displacements | ||
phonopy_atoms_to_ase_atoms(s) + atoms_to_add | ||
for s in phonon.supercells_with_displacements | ||
] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I need to return to this later to fully digest what's going on, but I'll trust you on it for the time being.
I don't really like the fact that |
I agree that it feels a bit wrong since the |
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## main #2073 +/- ##
=======================================
Coverage 99.06% 99.06%
=======================================
Files 81 81
Lines 3311 3325 +14
=======================================
+ Hits 3280 3294 +14
Misses 31 31 ☔ View full report in Codecov by Sentry. |
Last minute thought: "Fixing" atoms during phonons calculations is not a real thing and many people would simply call it blasphemy (to be fair |
@tomdemeyere: Thanks for the suggestion. I have been gravitating back towards having the My main hesitations right now are that |
"I have been gravitating back towards having the fix_indices keyword argument for this reason; the user really has to directly request it and know what they are getting themselves into." I think this might be better. When you suggested to use ASE's constraints I got hooked as it is definitely more consistent and "appealing", but after thinking about it a little bit more I would vote for keeping About Phonopy in Quacc: It might be more interesting to make simpler Phonopy functions instead of whole Users seeking something simpler might be better off being redirected to ASE's phonon/vibrations workflows (which are still powerfuls but arguably less complex). |
@tomdemeyere: Apologies again that this hasn't been as trivial to merge as I initially believed. I have opened up PR #2108 to demonstrate some areas for improvement. Namely, I want to draw your attention to There are still some problems with PR #2108, however:
Feel free to copy or pick things up from there. Unfortunately, I don't feel comfortable merging #2073 until this gets streamlined, especially for what amounts to a crude approximation in many cases (although one that is quite common).
Agreed. There is a constant, delicate balance between recipes doing too much and too little. At the end of the day, quacc has been intentionally designed to give advanced users all the tools they need to build their own custom workflows. Quacc isn't necessarily meant to codify every option available to users, and my vision is that end users can make their own scripts as they see fit. At the same time, when flexibility is available and makes sense to support, it is good to have. I don't think we need to try to overload the workflow to access every phonopy option available. |
I think I just got an idea that might avoid quite a lot of the current mess to fix atoms, Also avoiding get_phonopy returning a tuple. Let me have a look tomorrow. |
Closing in favor of #2110. |
Summary of Changes
I like the idea that you had, removing the parameter
fixed_atoms
and deal with constraints instead, some problems though: users might want to optimise these atoms before hand during therelax_job
, or sometimes not, how to deal with the constraints then? Also, this might be seen as a "hidden" feature?This MR also really removes the second
get_phonopy
which somehow evaded both of our attention during the last PR 😂Checklist
main
.Notes